Bubble and droplet motion in binary mixtures: 
Evaporation-condensation mechanism and Marangoni effect 
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Bubble and droplet motion in binary mixtures is studied in weak heat and diffusion fluxes and 
in gravity by solving the linearized hydrodynamic equations supplemented with appropriate surface 
boundary conditions. Without gravity, the velocity field is induced by evaporation and condensation 
at the interface and by the Marangoni effect due to a surface tension gradient. In pure fluids, the 
latter nearly vanishes since the interface temperature tends to the coexistence temperature T cx (p) 
even in heat flow. In binary mixtures, the velocity field can be much enhanced by the Marangoni 
effect above a crossover concentration c* inversely proportional to the radius R of the bubble or 
droplet. Here c* is usually very small for large R for non-azeotropic mixtures. The temperature and 
concentration deviations are also calculated. 

PACS numbers: 47.55.D-, 68.03.Fg, 64.70.F-, 44.35. +c 



I. INTRODUCTION 



On earth, bubble motion in liquid is caused by gravity. 
Buoyancy effects increase dramatically with increasing 
the droplet radius R. Due to the viscosity of liquid, it 
moves at a constant velocity v g estimated as [l|, 0, H| 



P-P 
V 



R 2 g, 



(1.1) 



where g is the gravity acceleration. Hereafter p and r\ 
{p' and rj') are the mass density and the shear viscosity 
outside (inside) the bubble. Another method of inducing 
bubble motion is to apply a heat flux Q. It is well-known 
that a surface tension variation on the surface gives rise 
to a Marangoni velocity field [l[, causing bubble motion 
to lower the surface free energy. Neglecting phase tran- 
sition, Young et al. [4| calculated it as 



Vy ~ —rRQ, 
1]A 



(1.2) 



where A is the thermal conductivity in liquid. The surface 
tension variation ^7 is assumed to be given by 



(1.3) 



where ST is the ambient temperature deviation, so vy ~ 
d'y/rj. Here 71 > for most fluids, but 71 < for some 
fluid mixtures. If a liquid is heated from a boundary 
at zero gravity, a suspended bubble is attracted to the 
warmer boundary for 71 > with a velocity of order 
— vy, until it is attached to the wall. In heat flux on 
earth, the gravity and Marangoni mechanisms can com- 
pete. We mention an experiment of applying heat flow 
from below to silicone oil containing air bubbles, where 
temperature gradients of order 10-10 2 K/cm balanced 
with the buoyancy and held the bubbles stationary 
Subsequent microgravity experiments on the Marangoni 
effect have been performed without phase change [5|, Q . 

However, first order phase transition between gas and 
liquid (evaporation and condensation) should take place 



on the bubble surface. This is particularly the case 
for pure (one-component) fluids, where the pressure p 
is nearly homogeneous outside the bubble for slow mo- 
tions and the interface temperature should then be close 
to the coexistence temperature T — T cx (j>) at given p 
even in heat flux. Thus, in pure fluids, the tempera- 
ture gradient should nearly vanish inside bubbles with- 
out Marangoni flow. Recently such temperature profiles 
have been calculated from linearized hydrodynamic equa- 
tions supplemented with appropriate surface conditions 
Q and numerically in the dynamic van der Waals theory 
8]. Balance of a heat flux due to latent heat convection 
and an applied heat flux Q gives the amplitude of the 
convective velocity inside the bubble as 



Q 



p'TAs ' 



(1.4) 



where As = s' — s is the entropy difference per unit 
mass. If a bubble (droplet) is suspended in liquid (gas) at 
zero gravity, it migrates toward a warmer (cooler) bound- 
ary with a velocity of order v c Q . In this evaporation- 
condensation mechanism, a bubble in liquid is attracted 
to a warmer boundary wall, which is consistent with ex- 
periments on pure fluids without gravity [9( . 

In this paper, we investigate bubble and droplet mo- 
tion in binary mixtures, where the Marangoni effect and 
the evaporation-condensation can be both important. In 
analyzing boiling experiments, Marek and Straub [1 0f | ar- 
gued that convection around a bubble should be domi- 
nantly caused by the Marangoni effect due to a very small 
amount of a noncondensable gas. If a surfactant is added 
as a solute, such a contamination effect should be even 
more enhanced [l], [ll[ . In our theory we shall see that 
the Marangoni velocity for dilute non-azeotropic binary 
mixtures is of order, 



VM 



7i 



k B nD rj 



cRQ, 



(1.5) 



where n is the number density, Do is the solute diffusion 
constant in liquid, and ks is the Boltzmann constant. 
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Balance of v c and vm yields a crossover concentration c* 
given by 



ai/R, 



(1.6) 



where a\ = (kBD r]/T)(n/ p' As\ji\) is usually a micro- 
scopic length far from the critical point. See sentences 
below Eq.(4.16) for a\ near the critical point. Unless 71 
is very small, c* is very small for large droplet radius 
R^> a\. For c>c* the hydrodynamic motion is mostly 
due to the Marangoni effect. 

Furthermore, there seems to have been no fundamen- 
tal argument on the coefficient 71 in Eq.(1.3) in binary 
mixtures in nonequilibrium. In this paper we assume the 
continuity of the temperature and the chemical poten- 
tials and neglect the pressure deviations at the interface. 
Then it follows 71 = ~(dj/dT) c ^ p , where the derivative 
is along the isobaric line on the coexistence surface. This 
71 tends to a well-defined limit in the dilute limit of bi- 
nary mixtures (c — > 0). Recently, it has been calculated 
for nonelectrolyte binary mixtures [l2| . Particularly near 
the solvent criticality, its mean-field expression at small 
solute concentration reads 



7i 



d-f dp c 1 
' dT dX K Kl 



(1.7) 



where 70 (T) is the surface tension of the pure solvent, 
Pc(X) is the critical pressure dependent on the solute 
molar fraction X on the critical line, and -Kxr is the so- 
called Krichevskii parameter (having the dimension of 
pressure) [lj, [Ij, [jja, [ly, [l7| • See Appendix C for discus- 
sions of K~Kr- While djo/dT < for pure fluids, the two 
parameters dp c /dX and can be both positive and 
negative, depending on the solute molecular size and the 
solute-solvent interaction. For example, if near-critical 
CO2 is a solvent, use of data in Ref.[17J gives the value 
of (dp c /dX)/ Kki in Eq.(1.7) for various solutes, which 
is 0.90 for neon but is —0.81 for pentanol [l2| . 

In their experiment, Vochten and Petre [H[ found 
that the surface tension between air and aqueous mix- 
tures containing high carbon alcohols more than ImM 
exhibits a minimum as a function of the temperature 
at constant pressure and molar fraction. In such fluids, 
71 < at temperatures higher than that giving the mini- 
mum. Inspired by their finding, consequences of negative 
71 have been discussed in two-phase hydrodynamics par- 
ticularly to develop heat pipes for utilization in space 
[H Ha [H H3, HI]. Remarkably, if the sign of ^ is 
changed, the direction of the Marangoni flow is reversed. 
As a result, if 71 < 0, bubbles are easily detached from 
the heater in boiling. This leads to a liquid inflow onto 
the heater suppressing its dryout, so fluid mixtures with 
71 < have been called self-rewetting fluids. On earth, 
Abe [2l[ observed a considerable decrease in the size of 
rising bubbles with addition of 1-butanol (6 wt%) to wa- 
ter (where 71 < 0). Adding 1-heptanol (0.1 wt%) to 
water, Savino et al. [23] observed bubble motion toward 
a cooler end in a horizontal glass tube. 



This paper will present linear analysis in the simplest 
case of a spherical bubble or droplet in binary mixtures in 
weak heat and diffusion fluxes and in gravity. In Section 
II, we will give linear hydrodynamic equations and sur- 
face boundary conditions including the Marangoni con- 
dition for the tangential stress [J [24j]. In Section III, 
we will solve the equations in steady states in the ax- 
isymmetric geometry. In Section IV, we will examine the 
consequences in dilute mixtures. Estimations near the 
critical point will also be presented. In Section V, the 
velocity field around a bubble or droplet will be displayed 
in various cases. 



II. BASIC EQUATIONS 
A. Spherical droplet 

We place a gas bubble in liquid with radius R in a non- 
electrolyte binary fluid mixture. We do not assume sur- 
face adsorption due to the amphiphilic interaction. The 
following results can be used also for the case of a liquid 
droplet in gas by exchanging "liquid" and "gas". Sup- 
pose an equilibrium state in the gravity-free condition 
(see Appendix A of Ref.[7]), where the temperature T 
and the chemical potentials of the two components /zi 
and /Z2 are homogeneous. The pressure p is a constant 
Po in the exterior r > R and is po + 2'y/R in the interior 
r < R ( from the Laplace law, where 7 is the surface ten- 
sion and r is the distance from the bubble center. The 
interior and exterior concentrations are determined from 
the thermodynamics of binary mixtures [l2| . 

We then apply weak heat and diffusion fluxes and a 
gravity acceleration g. They are all along the z axis taken 
to be in the upward vertical direction || e z . Hereafter e z 
denotes the unit vector along the z axis. The gradients 
of the temperature T and the mass fraction c are homo- 
geneous far from the bubble and are written as 



7 = (f 

dz 



, C=(|) • (2.1) 

r— 00 \ / r—00 



It is convenient to introduce the chemical potential dif- 
ference per unit mass as 



fj, = H2 - Mi, 
which has a gradient given by 



M \'± 

dz 



dTj pc 



dc ) 



(2.2) 



C, (2.3) 



pT 



Here the thermodynamic derivatives are taken in the 
outer phase in the isobaric condition. 

The hydrodynamic equations are linearized with re- 
spect to T, C, and g. The deviations are all proportional 
to one of these quantities. After a transient relaxation, 
the bubble moves at a constant velocity in the vertical 
z axis. We may then take the origin of the reference frame 
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at the bubble center and seek a steady axisymmetric so- 
lution of the hydrodynamic equations with appropriate 
boundary conditions. To linear order in T, C, and g, 
the bubble shape is spherical, as assumed in the previous 
theories 0, H, 0] • Deviation from sphericity occurs from 
second orders in these quantities. 

In the following calculation it is convenient to use the 
spherical coordinates (r, 9, tp) with the origin at the bub- 
ble center. Using the solid angles 9 and ip we define the 
three orthogonal unit vectors, 



where q is the dissipative heat current expressed as [2£ 



ei 
e 2 
e 3 



r r = (sin 9 cos p, sin 9 sin p, cos 9), 
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ei = (cos 9 cos ip, cos 9 sin <p, — sin 9) , 



= e\ x e 2 = (— sin ip, cos ip, 0). 



(2.4) 



The interface normal is along e\. The velocity field v 
around the bubble will be assumed to be orthogonal to 
e 3 or e 3 • v = 0. 



B. Hydrodynamic equations 

The mass densities of the two components are written 
as pi and p 2 . The total mass density p and the mass 
fraction c of the second component are defined by 



p = pi+p 2 , c = p 2 /p. 



(2.5) 



The mass fluxes of the two components are p\V — I and 
P2V + I, where v is the velocity field and J is the diffusive 
flux. The continuity equations for p and pi are [25[ 



d r—7 I \ 



-V • (pcv + I) 



(2.6) 
(2.7) 



We express J in terms of the isothermal mutual diffusion 
constant D and the thermal diffusion ration kx as [25j 



I = -pD(Vc + T _1 fc T VT). 
The momentum density pv obeys 

d «_> 

—pv + V • (pvv) = -Vp + V • cr - pge z 

at 



(2.8) 



(2.9) 



where p is the pressure and a = {cry} is the dissipative 
stress tensor expressed in terms of the shear and bulk 
viscosities 77 and C as 

(Tij = rj(Vi«j + Vw) + ((- 277/3)% V • v, (2.10) 

where Vj = d/dxi with Xi — x,y, and z. The last term 
in Eq.(2.8) represents the gravity acceleration with e z 
being the unit vector along the z axis. The (total) energy 
density ct = e + pv 2 /2 consisting of the internal energy 
density e and the kinetic energy pv 2 /2 is governed by 



_5 

<9t 



e T = -V • [(e T + p)w - a ■ v + q] - pgv z (2.11) 



-AVT + AI. 



(2.12) 



where A is the thermal conductivity in the absence of dif- 
fusion flux and A is a constant. Because of the symmetry 
of the Onsager coefficients (see Appendix A), there is a 
relation between fcy and A given by 



A 



t( — 



(2.13) 



From Eqs.(2.1) and (2.3) the heat flux and diffusion 
flux are written as q — ► — Qe z and J — » —Te z far from 
the bubble r ^> R, where 



Q = AT + .4J, 

J = pD(C + T- 1 k T T). 



(2.14) 
(2.15) 



The pressure gradient tends to —pge z , while the veloc- 
ity tends to —voe z because the bubble is at rest in our 
reference frame. In the present work we linearize the 
hydrodynamic equations for the velocity field v and the 
deviations Sp, ST, and 5c with respect to T, C, and g 
in steady states in the bulk region r 7^ R. Here we may 
set d(- ■ ■) / dt = neglecting the time-dependence. From 
Eq.(2.6) it follows the incompressibility condition, 



V • v = 0. 



(2.16) 



Then the bulk viscosity £ does not enter in our calcu- 
lations. This incompressibility condition holds even in 
compressible fluids in our linear theory. Notice that we 
neglect the term v ■ Vp in Eq.(2.6) since it is of order Tg 
in the gravity-induced density stratification. In the bulk 
region r ^ R, Eq.(2.7) and (2.11) yield 

V 2 5T = 0, V 2 (5c = 0, (2.17) 

The momentum equation becomes 

- V5p + r]V 2 v - pge z = 0. (2.18) 

The pressure deviation is defined by Sp = p—po for r > R 
and Sp — p — po — 270 /R for r < R. For r ^ R, taking 
the divergence of Eq.(2.18) yields 



V 2 Sp = 0. 



(2.19) 



The term —(dp/dz)g is of the second order for r ^ R 
and is negligible in our approximation. 



C. Interface boundary conditions 

Next we consider the boundary conditions at the sur- 
face r = R. In this subsection all the quantities are those 
at r = R ± 0. In the following equations the quantities 
at r = R — (immediately inside the bubble) are primed 
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as v', ST', Sp' ■ ■ ■, while those at r = R + (immedi- 
ately outside the bubble) are unprimed. Hereafter, for 
any physical quantity A, the symbol, 

[A]=A-A', (2.20) 

denotes the discontinuity of A at the surface. For exam- 
ple, the entropy difference (per unit mass) and the mass 
concentration are written as 

[s]=s-s', [c] = c-c'. (2.21) 

The Gibbs-Duhem relation yields [s]/[c] = — (dp/dT) CXtP 
in terms of the chemical potential difference p. Hereafter 
{&■■■/&■■ -)cx,p is the derivative taken along the isobaric 
line on the coexistence surface. 

The mass conservation at the surface yields e.\ ■ [pv] = 
0. It is convenient to introduce the mass flux J through 
the interface by 

J = pei ■ v = p'ex ■ v', (2.22) 

which arises from conversion between gas and liquid. We 
assume the continuity of the tangential velocity, 

e 2 • [v] = e 2 • v - e 2 • v' = 0. (2.23) 

The stress balance at the interface yields 

[8p-e 1 -V-e 1 } + ^7 = 0, (2.24) 
[e a • a • e x ] + e 2 • V<5 7 = 0, (2.25) 

in the normal and tangential directions, respectively. 
Here £7 = 7 — 70 is the surface tension deviation and 
e 2 • V = Rr x d/d6. From Eq.(2.25) the tangential gra- 
dient of the surface tension is equal to the discontinuity 
of the t ang ential stress, which gives rise to a Marangoni 
flow 0, . 

As in the pure fluid case 0, we assume that the de- 
viations of the temperature and the chemical potentials 
are continuous at the interface, 

ST = ST', (2.26) 
Sp± = Sp[ <5/j 2 = 5p 2 - (2.27) 

See Appendix B of our previous work Q and the item 
(i) in the summary section of this work for discus- 
sions on the validity of these assumptions. The Gibbs- 
Duhem relation for infinitesimal deviations is written as 
(1 — c)5p\ + cS/i2 = —sST + p~ x 5p, which holds in the 
liquid and gas regions close to the surface. Further use 
of Eqs.(2.26) and (2.27) yields 

[c]Sp + [s]ST = -Sp - —Sp', (2.28) 
p p' 

where Sp = Sp2 — Sp\. 

From Eq.(2.7) the mass conservation of the second 
component at the surface gives 

[c]J + ei ■ [/] = 0. (2.29) 



From Eq.(2.11) the energy conservation at the interface 
gives 

([ S ]T+[c]p)J + ei -[q}=0, (2.30) 

where use has been made of the thermodynamic relation 
e + p = pipi + P2P2 + psT = (pi + cp + sT)p. From 
Eqs.(2.29) and (2.30) J may be removed to give 

e x -[q- aTI] = 0, (2.31) 

where the coefficient a is defined by 

a = p/T+ [s]/[c]. (2.32) 

Note that a is a constant continuous across the interface. 
The flux q — aTI is continuous along the normal direction 
across the interface. 



D. Pressure and surface tension deviations 

We shall see that the pressure deviations in the two 
phases are negligibly small for large R in Eq.(2.28). This 
yields the following relation, 

which plays a key role in the following calculations. It 
may be justified if the hydrodynamic deviations are ex- 
panded in powers of the inverse bubble radius R . To 
leading order in R^ 1 , the right hand side of Eq.(2.28) is 
negligible, resulting in Eq.(2.33). 

In equilibrium, the surface tension of binary mixtures 
is defined on the coexistence surface in the space of three 
field variables such as T, p, and p. Here the curvature 
effect is neglected for large R. Then the surface tension 
is a function of T and p if p is taken to be the coexis- 
tence pressure p cx (T,p). In Eq.(2.33) the temperature 
and chemical potential deviations near the interface are 
still on the coexistence surface in the isobaric condition. 
Thus the surface tension deviation is written as [27] ] 

<5 7 = a\ST + a 2 <5/i 

= -7i<5T, (2.34) 

where a\ and a 2 are the expansion coefficients. As dis- 
cussed in Section 1, the coefficient 71 is written as 

III. AXISYMMETRIC SOLUTION 

A. Velocity and pressure 

In our problem, the fluid flow is axisymmetric with 
respect to the z axis. The velocity field v(r) and the 
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pressure deviation 8p(r) are expressed in the same forms 
as in the one-component fluid [7, 26]. That is, in terms 
of two functions Q(r) and H(r), v(r) is written as 

dH H ~\ z 1 - . , 

v I + Qr)^r + -He z . (3.1) 

ar r J r z r 

which satisfies Eq.(2.16). Outside the bubble we have 

Q = Qx^, H=^Q 1 + H 1 ^~v D r (r>R). (3.2) 

where «d is the bubble velocity in the original reference 
frame. Inside the bubble we have 

Q = Q' 2 ^, H = -^Q' 2 r 3 + H 2 r (r < R). (3.3) 

Then Eq.(2.14) is satisfied. In particular, Eq.(2.16) yields 

Qi=gR 2 (p-p')/3r,. (3.4) 

The coefficients Qx, Hi, Q' 2 , and H 2 have the dimen- 
sion of velocity. From Eqs.(2.18) and (2.19) the pressure 
deviation is determined as 



Sp(r) = vQi^r ~ 9P Z 



(r > R) 



= -^Q' 2 —~gp'z (r<R). (3.5) 

The mass flux through the interface and the surface 
tensition deviation are angle-dependent as 



J = JiCos6>, Sj — YxCOsO, 



(3.6) 



where cos# = z/r and Ji and Yi are constants. From 
the boundary conditions Eq.(2.22)-(2.25), Hi, Vr>, H 2 , 
and Ti may be expressed in terms of the three quantities 
Ji, Q' 2 , and Q\ as 

Hi = \{^--)Ji-^Q' 2 + \Qu (3.7) 
3 p p 15 6 

«d = -\{^ + -)Ji + ^Q'2 + \Qi, (3-8) 
3 p' p 15 3 



H 2 — —Ji + -Q' 2 , 
p' 5 



(3.9) 



Ti = (^ - ~)vJi - &T) + \rf)Q' 2 + r,Q x . (3.10) 
p p 5 5 

Since Qi is determined as in Eq.(3.4), the two quantities 
J i and Yi remain unknown. 

From Eqs.(2.18) and (2.19) the deviation of the chem- 
ical potential difference dp, also satisfies V 2 Sp = for 
r ^ R. Then ST and Sp are written as 

5T(r) = (T' ~T)^z + Tz (r > R) 

= T'z {r<R), (3.11) 

R 3 

5p(r) = (M - M)—z + Mz (r > R) 



where T and M. are the gradients far from the bubble 
defined in Eqs.(2.1) and (2.3), while T and M! are those 
within it. At r — R, ST and Sp are continuous from 
Eqs.(2.26) and (2.27). For r < R we have 

ST = T'z, Sp = M'z. (3.13) 

In terms of T and M! the relation Eq.(2.28) is rewritten 
as 



\c]M' 



]T' = 



-Qi + 2 4Q' 2 
p pf 



Furthermore, we may derive two equations for T' and 
M' from Eqs.(2.29) and (2.30). Together with Eqs.(3.10) 
and (3.14), we have four equations, which constitute a 
closed set determining J\, Q' 2 , T' , and M! . Here, without 
assuming Eq.(2.33), we may assume the general relation 
67 = a x 8T + a 2 Sp in the first line of Eq.(2.34). 



B. Case of g ^ and Q = X = 

Here we consider the gravity-induced solution with 
g =^ and Q = T — 0, where all the coefficients in 
Eqs.(3.7)-(3.10) are proportional to Q x in Eq.(3.4). If 
the pressure deviations in Eq.(2.28) are neglected, there 
arise no deviations of the temperature, the chemical po- 
tential, and the surface tension, T = T' = J x = Yi = 0. 
in Eqs.(3.7)-(3.10). Then Eq.(3.10) gives 



Q' 2 = 5r,Qi/(2r, + 3 V '), 



(3.15) 



while Eq.(3.8) gives vd = Va, where v g is the well-known 
gravity-induced velocity 0,13, E| . 



2(?7 + r ? 0(p-pQ 
377(277 + 377') 9 ' 



(3.16) 



In the lowest order of R 1 the right hand side of Eq. (2.28) 
at r — R becomes 



Sp Sp' 
p p' 



107777' 



n 



p p'(277 + 377') 



R 2 



(3.17) 



M'z 



(r<R), (3.12) 



where the gravity terms cancel to vanish. If divided by 
Qx, the above quantity is of order If we assume 

the linear relation $7 = axST + a 2 Sp, the above relation 
Eq.(3.17) leads to ST oc Spot Qx/R, T' oc M' oc Qx/R 2 , 
Yx oc Qx/R, and Jx oc Qx/R 2 . Thus, in the presence of 
weak gravity only, a large droplet or a large bubble moves 
with the velocity v g in Eq.(3.16), where first-order phase 
transition and temperature inhomogeneities are negligi- 
ble. 



C. Case of Q / 0, I + 0, and g = 

We seek the solution in the presence of Q and X in 
the gravity- free condition g = 0. Remarkably, Yx oc R 
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in binary mixtures with phase change, while Ti oc 
in one-component fluids. We clarify the relationship of 
our theory and the previous theories: (i) To obtain the 
solution without phase change [J], we set J\ = and 
ri = -7ii?T oc R in Eqs.(3.7)-(3.10) and require the 
energy conservation relation e\ ■ [q] = from Eq.(2.30). 
(ii) To obtain the solution for one-component fluids with 
phase change [7[, we neglect Ti and retain J\ in Eq.(3.10). 

We use the relation Eq.(2.33) or neglect the right hand 
side of Eq.(3.14) to obtain 



M' = -T'[s]/[c] 



(3.18) 



The above relation will be justified self-consistently at 
the end of this subsection. From Eqs.(2.34) and (3.13) 
Ti in Eqs.(3.6) and (3.10) is expressed as 

r : = -7X-RT'. (3.19) 

As will be shown in Appendix B, T can be written as 

3 



T = 



2A e + X' e 



7 {Q~aTl), 



(3.20) 



where Q and X are defined by Eqs.(2.14) and (2.15). We 
introduce the effective thermal conductivity by 



A e = A + 



zf 



(3.21) 



The A e and \' e in Eq.(3.20) are the values of A e at r = 
R + ±0. In Eq.(3.21) we define Z by 



Z 
f 



dc_ 

df 



dc 

dfx ) 



Tp I 



dc 
df 



(3.22) 



On the other hand, the mass flux through the interface 
J is calculated from Eq.(2.29) or Eq.(2.30). Then J x in 
Eq.(3.6) is expressed as 



Ji 



< (2X a + X' a )l-(2B + B / )Q 
[c](2A e +A'J 



(3.23) 



where A a , X' a , B, and B' are the values of A a and B at 
r = R + ±0 with 



A a = A e + aTB = A + AB, 
B = pD(k T + Z)/T. 



From Eq.(3.10) we obtain 
5 



Q' 2 = 



2ri + 3rj' 



7l Rr + (1 - -)vJi 
P P 



(3.24) 
(3.25) 



(3.26) 



Substitution of the above relation into Eq.(3.8) yields the 
drift velocity composed of two parts, 



t'D 



(3.27) 



where arises from the Marangoni effect and v c D from 
the evaporation-condensation. They are written as 



M 



I'D 



2l 



3(2^ + 3r}') 



(3.28) 



- 2 Jl + &i)J>- ("») 



2rj + 3rj' \p 2p p 1 



In terms of these characteristic velocities, the velocity 
field v(r) in the reference frame moving with the bubble 
or the droplet is expressed as 



M 

v = v D a M - 



V(p' - P)v c D 



i 1P > + rf(p + p>/2) 



L—a c -(v%+v c D )e z , (3.30) 



where ajj = aif(r) and a c = a c (r) are the following 
space-dependent dimensionless vectors, 



2r3 6 



3R j 



zr (r > R) 



2r) 
zr 



2r 5 

% Bz+ 2l? Zr (r<i?) ' (3 - 31) 



3R 3 



<1 



2r 
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2rj R 2 



(r > R) 



e z (r < R), (3.32) 



We may now show that the right hand side of Eq.(3.14) 
is surely negligible for large R. Use of Eq.(3.26) gives 



Sp Sp' _ 10r/ 
p p' 2i] + 3r/ 



7i T' 2 2 yJ, 
p'R V p'p'R 2 



(3.33) 



We compare the term (cx i? _1 T') on the right hand side 
of Eq.(3.33) and the term [s]T ! on the left hand side of 
Eq.(3.14). The former is much smaller than the latter 
for R 3> |7i/p'[s]|, where the right hand side is micro- 
scopic. The term (cx R~ 2 Ji) due to the evaporation- 
condensation in Eq.(3.33) is also negligible, as already 
verified in our previous paper [7[. 

The concentration deviation Sc(r) can also be ex- 
pressed in the same form as in Eqs.(3.11) and (3.12), 
where the gradient C far from the droplet satisfies 
Eq.(2.3), Its values at r = R ± are written as 5c± and 
are expressed as 



5c, 



T 



T'z, Sc- 



~T 



T'z, 



(3.34) 



where Z and Z 1 are the values of Z in Eq.(3.22) at r — 
R±0. If Z / Z' , 5c is discontinuous at r = R. 

We note that the Marangoni effect vanishes in 
azeotropic mixtures 26], where the two phases have the 
same composition or [c] =0. There is also no difference 
in the molar fractions of the two phases (see the sentence 
below Eq.(C3)). Special analysis is thus needed when we 
treat nearly azeotropic mixtures. For example, in H 2 0- 
D2O mixtures, the relative composition change [c]/c is 
only 0.5% of the relative density change [n]/n near the 
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critical line [l7|. Thus let us consider the limit [c] — > 
in the equations in this subsection. Then Z ~ [c] _1 
and A e ~ [c]" 2 from Eqs.(3.20) and (3.21), leading to 
r x ~ 7l [c] 2 , V ~ [c] 2 (Q - aTI), and 



7i[c] 2 i?(Q-«rc). 



(3.35) 



from Eqs.(3.19), (3.20), and (3.28). Here 71 - [c]" 1 [T^ 
and 71 [c] 2 ~ [c], so the Marangoni flow is of order [c]. 



where Z?o = lim c ^o -D is the diffusion constant of a single 
solute molecule and n — p/mi is the solvent number den- 
sity. The thermal conductivity A in Eq.(3.21) is smaller 
than the right hand side of Eq.(4.7), as will be discussed 
in Appendix C. Here, in calculating the flux e.\ -(q — aTI) 
in Eq.(2.31), we have picked up the contribution from the 
solute diffusion to obtain A e in Eq.(4.7). 

Furthermore, A a in Eq.(3.24) and B in Eq.(3.25) are 
finite as c — > 0. From Eq.(3.20) we thus find T oc c as 



IV. DILUTE MIXTURES IN GAS-LIQUID 
COEXISTENCE 



Let the second component be a dilute solute. Under 
the condition Eq.(A8) we set 



(dc/dfi)rp — m 2 X/mxT, 



(4.1) 
(4.2) 



where X is the molar fraction and mi and m2 are the 
molecular masses. The second relation (4.2) does not 
hold very close to the solvent criticality even for small 
X (see Eq.(C7)). Hereafter the Boltzmann constant will 
be set equal to unity. In the literature [12, [l3|, 0, HE 
liH \v\ , Henry's law is expressed in terms of the solute 
molar fractions. That is, in equilibrium, the solute molar 
fraction in gas X g and that in liquid Xi are related by 
the partition coefficient [12j, 



K = Xg/Xe, 



(4.3) 



which depends on T along the solvent coexistence line 
p = Pcx(T). In dilute mixtures [c]/c = [X]/X holds and 
[c]/c is independent of c as 

[c]/c = 1 — K, (gas bubble) 

= 1-/CT 1 (liquid droplet). (4.4) 



A. Expressions as c — > 

To simplify the notation we introduce the following 
dimensionless parameter, 



W 



[X] _ [c] 
X[a] ctoi[s]' 



(4.5) 



where a = m\s is the entropy per solvent particle. This 
parameter tends to a well-defined limit in the dilute limit, 
vanishes for azeotropic mixtures, and becomes propor- 
tional to the Krichevskii parameter Kk v near the solvent 
criticality as in Eq.(C6) in Appendix C. In Eq.(3.22) Z 
behaves in terms of W as 



ni2c[s] 



m 2 



[c] miW 
From Eq.(3.21) A e is inversely proportional to c as 



A,. 



m 2 c 



W 2 A' 



(4.6) 



(4.7) 



T'^A M c Q-T^I 



(4.8) 



Here we have set a = [s]/[c] from Eq.(2.32) assuming 
that X is of order c. The coefficient Am is a positive 
constant independent of c as c — > and is defined by 



A 



M 



3W 2 mi/m2 
2nD + X'u'D'q/X 1 



(4.9) 



which is expressed in terms of the molar fractions X and 
X' outside and inside the domain. The mass flux through 
the interface J in Eq.(2.22) tends to that of the pure fluid 
as c — > 0. That is, Ji in Eq.(3.6) becomes 



Ji = (3/WT)Q. 



(4.10) 



Now the two drift velocities v|f in Eq.(3.28) and v c D in 
Eq.(3.29) are written as 



M 



2AM11 cR(Q- [ ^Tl), (4.11) 



3(2?7 + 3r?' 

— 6 / 7] 

2r) + 3r]' \p 



2 P P ')[s\r 



(4.12) 



From Eq.(3.34) the deviation of the mass fraction 5c 
at r = R ± are expressed as 



c 



Sc. 
c 



^ ± = ^ 1 _ = m 2 [s\ , 



(4.13) 



Thus 5c/ c is continuous and [5c] — — rri2[s\5T at r = R 
to leading order in R . In our linear theory we re- 
quire \5c\ <C c, which becomes |T'|i? <C |[c]/m 2 [s]|. From 
Eq.(4.8) this inequality is satisfied for small Q and X even 
as c — > 0. 

Let us consider situations in which the diffusion flux 
X is negligible in in Eq.(4.11). Then v 1 ^ and v c D 
have the same sign for a bubble with 7l > and for 
a liquid droplet with 71 < 0, while they have different 
signs for a bubble with 71 < and for a liquid droplet 
with 71 > 0. In accord with the experiments [2ll. I22I [23]] . 
bubbles can move towards cooler regions with increasing 
the concentration of a solute in the case 71 < 0. 



B. Dilute mixtures near the solvent criticality 

The above expressions can be used even in the vicin- 
ity of the solvent criticality under the condition X <C 
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n/W 2 C p in Eq.(C8), where C p is the isobaric specific 
heat per unit volume. In the near-critical case, the re- 
duced temperature, 

e = 1 - T/T c , (4.14) 

is small (say, less than 10 -3 ) and the differences between 
the two phases tend to vanish, so p' = p and rf = -q. In 
Eq.(4.9) we have 

A M = \\> 2 m 1 /m 2 nD Q . (4.15) 

Defining the solute hydrodynamic radius do using the 
Stokes formula Dq = T /Qitrjao, we obtain 

v ° = i^ w2liaoRX { Q - [ i] TI )> (4 - 16) 

v c D = -3Q/nT[a], (4.17) 

where [a] — mi [s] ~ e 13 . 

With increasing X from zero, crossover occurs from the 
pure-fluid behavior to the mixture-behavior for X > X* , 
where the crossover molar fraction is 

X* = {e p \ 11 \W 2 a y 1 R- 1 . (4.18) 

This is equivalent to Eq.(1.6) near the solvent criticality 
withe* = m 2 X*/ mi ~ X* if |W| ~ 1 and m 2 ~ m x . The 
right hand side is of order e 1_!y_ ^^/i? for |W| ~ 1 and 

7i djo/dT with v = 0.625 and (3 = 0.33, £(~ e"") 

being the correlation length (~the interface thickness) . 

In near-critical pure fluids in the gravity-free condi- 
tion, a bubble in liquid was observed to be attracted to 
a warmer boundary 0] . With addition of a small amout 
of various solutes, it is then of great interest whether a 
bubble is more attracted to or eventually repelled from a 
warmer boundary. Here the crossover concentration X* 
should be measured to confirm the theoretical expression 
(4.18). 



V. VELOCITY PROFILES FOR Q ^ 0, 1 ^ 0, 

AND g = 

Young et al. [4] calculated the velocity field without 
phase change due to the Marangoni effect for 71 > 
(which is given by Eq.(3.30) with v c D = 0). See its dipolar 
profile in their paper and in our previous paper On 
the other hand, our previous paper Q has presented some 
examples of the velocity field due to first-order phase 
change for pure fluids (which is given by Eq.(3.30) with 
= 0). Here we display the velocity fields realized in 
dilute mixtures with increasing c or in nearly azeotropic 
mixtures. 

Without gravity, we show the velocity field (v x ,v z ) in 
the x-z plane for a bubble with p' j p — Tj'/rj — 0.5 in Fig.l 
and for a liquid droplet with p 1 / p = rj /rj = 2 in Fig. 2. It 
is written in the reference frame moving with a bubble 



or a droplet in the case v c D < 0. Here the liquid density 
is twice larger than the gas density, which is realized for 
T/T c = 0.97 in the van der Waals theory of pure fluids. 
In the left panels, we set = for pure fluids or for 
71 = in binary mixtures. In the middle panels, we set 
v d = v< b' wnere the two mechanisms equally contribute 
to the drift velocity. In the right panels, we set = 
—Vjj, where the drift velocity Vu in Eq.(3.27) vanishes 
and the bubble or the droplet is at rest. 



VI. SUMMARY AND REMARKS 

We have examined the competition of the evaporation- 
condensation effect and the Marangoni effect in the mo- 
tion of a bubble or a droplet in weak heat and diffu- 
sion fluxes in binary mixtures. We have treated the 
simplest case of steady states with a constant drift ve- 
locity in the axisymmetric geometry, though the non- 
linear terms in the hydrodynamic equations cannot be 
neglected in practical applications using large bubbles. 
In non-azeotropic binary mixtures, the crossover occurs 
from the evaporation-condensation mechanism to the 
Marangoni mechanism at a very small solute concentra- 
tion. In our theory, the coefficient 71 in Eq.(1.3) controls 
the strength of the Marangoni flow in heat flux, which can 
be both positive and negative depending on the solvent 
and solute species [ijj]. The Marangoni flow is induced 
in opposite directions in the normal case 71 > and in 
the anomalous case 71 < 0, though the case 71 > has 
mostly been studied. 

Some further remarks are given below. 

(i) We have assumed the continuity of the temperature 
and the chemical potentials and neglected the pressure 
deviations at the interface to obtain the key relation 
Eq.(2.33) or Eq.(3.18). It means that the interface 
stays on the coexistence surface p — p cx (T, /x) even in 
nonequilibrium. This is justified for large R. However, 
if the gas phase is very dilute far below the critical 
temperature, the surface dissipation mainly occurs in the 
gas phase side within a distance of the mean-free path 
inversely proportional to the gas density [H, [29|, [30| . 
There can then be an apparent temperature jump at the 
interface. 

(ii) The behavior of the coefficient 71 in Eq.(1.3) or in 
Eq.(2.35) is highly nontrivial. Theoretically, it has been 
examined only for dilute mixtures [l2l ]. Its behavior is 
also of interest for binary mixtures near a lower critical 
solution temperature (LCST). In a phase-separated 
mixture of butoxyethanol-water near its LCST, Braun 
et al. [3l[ applied heat-pulses to water-rich droplets to 
observe their motion from a high-temperature region 
to a low-temperature region. This motion was due to 
the Marangoni effect because the motion was in the 
direction of decreasing the surface tension. 

(iii) Boiling on a heated substrate has been of great 
interest both on earth and in space [32]. The effect of 
a noncondensable gas should be studied in future. In 
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FIG. 1: Gravity- free velocity field around a bubble for 77 /r) 
(left), v% = v c D (middle), and v% = -v° D (right). 




p / p = 0.5 in the x-z plane in Eqs.(3.30)-(3.32) 



where Vjj 
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FIG 



2: Gravity-free velocity field around a liquid droplet for rj jr\ — p / p = 2 in the x-z plane in Eqs.(3.30)-(3.32), where 



M = (left), 



(middle), and v£ = -v c D (right). 



accord with this paper, Marek and Starub [10( claimed 
that the temperature gradients along the bubble inter- 
face inducing a Marangoni flow are caused by saturation 
pressure gradients due to a nonuniform accumulation 
of a noncondensable gas along the interface. Such flow 
serves to suppress detachment of bubbles for 71 > 0, but 
should accelerate it for 71 < plj . 

(iv) In gravity-free conditions, a spherical bubble or 
droplet can be suspended in liquid or gas in equilib- 
rium. It is of great interest how the velocity field 
and temperature evolve after application of heat flux 
from a boundary. The piston effect comes into play on 
acoustic time scales [26|, |33 |. On longer time scales, 
a small amount of a solute should drastically change 
the hydrodynamic behavior inducing a Marangoni flow. 
However, since we have treated only steady states, it 
remains unclear how the concentration changes in time 
along the interface. 

(v) Thermocapillary hydrodynamics has been puzzling 



near the critical point H^, [33[ , where the singularities 
of the thermodynamic and dynamical properties largely 
influence the dynamics. The condition (4.2) or the 
condition (C8) does not hold sufficiently close to the 
critical point, where the results in Section 4 cannot be 
used. Thus, with addition of a solute, two-phase hydro- 
dynamics poses a new problem of critical dynamics, 
(vi) Finally, we should stress that surfactant molecules 
absorbed at interfaces give rise to Marangoni flow [l], [ll[ , 
though this effect is beyond the scope of this paper. 
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Appendix A: Onsager relations 

In the linear response theory, the dissipative heat and 
diffusion fluxes, q and J, in binary mixtures are expressed 
in terms of the Onsager kinetic coefficients L^ [25[ ■ The 
thermodynamic forces are VT and V(/i/T) as 



q = -L n VT-TL 12 v|, 
I = -L 12 ^VT-L 22 \7^, 



(Al) 
(A2) 



where p — p 2 — pi with pi being the chemical potentials 
per unit mass. We then consider the dynamic equation 
for the entropy density S = ps per unit volume. The 
thermodynamic relation TdS — de — p\dp\ — p 2 dp 2 and 
the hydrodynamic equations yield 



d ( 1 



p_ 

T J 



1 

= T 6 ' 



(A3) 



where e is the heat production rate per unit volume ex- 
pressed as 



e = J> 



13 dxi 



^-TI-V£. (A4) 



The first term arises from the viscous damping and the 
last two terms from the heat conduction and diffusion. As 
is well-known, a is nonnegative-definite if the coefficients 
Lij constitute a symmetric positive-definite 2x2 matrix. 
It is well-known that are expressed in terms of the 
the time-integral of the appropriate flux time-correlations 
(Green-Kubo formulas) [26J. 

It is convenient to express q as in Eq.(2.12) and I as 
in Eq.(2.8). Then A, A, D, and &t are expressed as 



A — in — L\ 2 / L 22 , 
A = TLi 2 /L 22 , 

dp\ 



pTD = L 22 
pDk T = Li2 



9c J pT 



T 



dp 
df 



(A5) 
(A6) 

(A7) 
(A8) 



cp 



If Li 2 and L 22 are removed from Eqs.(A6)-(A8), kx and 
A are related as in Eq.(2.13). 

Appendix B: Calculations of T and J\ 

To derive Eq.(3.20) we calculate the heat and diffusion 
fluxes at the interface in the normal direction substitut- 
ing Eqs.(3.11) and (3.12) into Eqs.(Al) and (A2) and 
setting M' = -T[s]/[c] from Eq.(3.14) for large R. The 
unprimed quantities are the values at r = R + 0, while 
the primed ones are those at r = R — 0. 

From Eq.(3.11) the gradient e.\ ■ VST normal to the 
interface is (3T - 2T') cos 9 for r = R + and to T cos 6 
for r — R — 0. From Eq.(3.12) the gradient e\ ■ V5p 



is obtained by replacement of T' and T by M! and M. 
Then use of Eqs.(3.11)-(3.14) gives 

e x q = (2A Q T'-3Q)z (r = R + 0), 

= -X' a T'z (r = R-0), (Bl) 

ei / = (2BT' - 3J)z (r = R+ 0), 

= —B'T'z (r = R-0) (B2) 

where z = z/r — cos9. We define Q and X in Eqs.(2.14) 
and (2.15) and introduce 

A a = L n -aL 12 = X + L 12 (A/T -a), (B3) 
B = {L a - aL 22 )/T = L 22 {A/T - a)/T. (B4) 

The flux q — Tal along e\ in Eq.(2.31) may then be 
calculated at r = R ± 0. Its continuity at r = R gives 

(2A e + X' e )T' - 3(Q - Tal) = 0. (B5) 

Bere A e is the effective thermal conductivity defined by 

A e = A a - aBT 

= \ + L 22 (A/T-a) 2 , 

where use of Eqs.(2.13) and (2.32) yields 

.4 



(B6) 



7^ — Ci 

T 



dp \ kr _( dp}\ _H 
dc) Tp T [dT y 



Cp 



(B7) 



From (dc/p) T p(dp/dT) cp = -(dc/dT)„ p and Eq.(A7) 
we obtain Z in Eq.(3.22) and A e in Eq.(3.21). From 
Eqs.(2.29) and (3.6) J\ in the mass flux through the in- 
terface is determined by 

Ji - [3fi - (2A„ + X' a )T']/T[c]a 

= [31- (2B + B')T']/[c], (B8) 

which yields Eq.(3.23) with the aid of Eq.(3.20). 

Appendix C: Mass and molar fractions and 
relations in near-critical dilute mixtures 

In hydrodynamic theory the mass fraction c and the 
chemical potential difference p per unit mass are usually 
used, but in thermodynamics it is convenient to use the 
molar fraction X and the chemical potential difference A 
per particle. In terms of the molecular masses mi and 
m 2 they are related by [26j 



c = m 2 X/[mi(l -X) +m 2 X], 
A = m 2 p 2 - m\pi. 



(CI) 
(C2) 



At any concentration it generally holds the relation, 



I^J =m 2 1 m 2 2 (n/p) s (^- 



Op 



TP 



OA 



(C3) 



TP 
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We also note the relation [c] = \X\[m\cd jm^XX'). The 
azeotropy condition [X] = is also given by [c] = 0. 

First we give thermodynamic relations in dilute mix- 
tures with c <C 1. In the text of this paper we assume 
Eqs.(4.1) and (4.2). To describe the critical behavior, it is 
convenient to introduce the Krichcvskii parameter Kkt, 
which is the dilute limit of the thermodynamic derivative 
(dp/dX) n T at fixed n = ni+n2 and T at the solvent crit- 
icaUty pi [ii [ii [Mini- 

It is related to the derivatives 
of the critical pressure p c and the critical temperature T c 
with respect to X along the critical line as 



dp c 
dX 



l dX' 



(C4) 



where p' cx = dp cyi /dT is the temperature-derivative of 
the coexistence pressure p cx (T) of the pure fluid. The 
Clausius-Clapeyron relation p' cx — [<j]/[v] holds for the 
pure fluid, where v = 1/n is the inverse density and a is 
the entropy per particle. In near-critical two-phase coex- 
istence, the mass fraction difference [X] and the volume 
difference [v] are related by [l2j 



[X]/[v] = {K Kl /T c )X c 



(C5) 



where X c = (X g + Xt)/2 is the critical value of X [T3|. 
Thus the parameter W in Eq.(4.5) is expressed near the 
criticality as 



W = K Kl [v]/T c [a] = K Kl /T c p' cy 



(C6) 



On the other hand, the the thermodynamic derivative 
(dX/dA) Tp behaves as [H 



dX \ ^ X 
dAJ ~ ~ 



pT 



X^ 



T nT 2 



K c1 .K T a- 



(C7) 



where Kta = {dn/dp)TA/n is the isothermal compress- 
ibility growing strongly near the criticality. On the right 
hand side of Eq.(C7), the first term is the dilute limit, 
while the second term is the singular contribution stem- 
ming from the solute-solvent interaction and can be im- 
portant very close to the criticality. In this paper we 
neglect the second term in Eq.(C7), which is allowable 
under the condition, 



nT 



K 2 K T W 2 C P 



(C8) 



where Kt = (dn/dp)T/n is the isothermal compressibil- 
ity and C p = pT(ds/dT) p = T(p' cx ) 2 K T is the isobaric 
heat capacity per unit volume of the pure fluid [2(| • How- 
ever, the reverse relation X > n/W 2 C p eventually holds 
sufficiently close to the criticality. 

Next we consider the dynamic properties of dilute mix- 
tures. The Onsager coefficients Lyi and L22 are propor- 
tional to c, while Ln tends to the thermal conductivity 
of the pure fluid. Thus, as c — > 0, A in Eq.(A6) and D in 
Eq.(A7) tend to well-defined limits, while h? in Eq.(A8) 
is proportional to c as 

k T = k* T c, (C9) 
with fcy being a constant independent of c. The singular 
part of L22 (proportional to X 2 ) is ne gligi ble compared 
to the background part under Eq.(C8) [2y|, so that 



m 2 pcD , 



(CIO) 



where Dq is the diffusion constant of a single solute par- 
ticle in the dilute limit. The thermal conductivity A in 
Eq.(A5) behaves in a dilute binary mixture near the sol- 
vent criticality as [2(| 



A = A S A B /(A S + A B ). 



(Cll) 



Here A s is the thermal conductivity of the pure fluid 
growing strongly near the solvent criticality and As = 
Xbo/X is the critical value inversely proportional to c 
with Xbo being a constant. The mode-coupling the- 
ory of critical dynamics predicted the growing behavior 
X s = TCp/Girrjt; ~ £, where £ is the correlation length 
and Cp{^ Kt) is the isobaric heat capacity per unit 
volume. We then recognize that the effective thermal 
conductivity A e in Eq.(3.21) may be approximated as 
Eq.(4.6) under Eq.(C8). Using Eq.(C8) and the hydro- 
dynamic radius ao in the Stokes formula Do = T/6irnaQ, 
we can make the following estimation, 



A 



> 



A, 



nDp 
W 2 XK 



n£/a 
C P W 2 X 



> 1. 



(C12) 
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